rm(list = ls())
library(tidyr)
library(readr)
library(dplyr)
library(readxl)
library(tidyr)
library(arm)
library(lme4)
library(MASS)
library(texreg)
load("Data/paperdata.RData") #1993-2006
# paperdata <- paperdata %>% dplyr::filter(gwno_loc !=710)
# save(paperdata, file = "Data/paperdata.RData")
source("Code/coefplot.R")

paperdata <- paperdata  %>% 
  dplyr::arrange(conflict_id,dyad_id, year) %>% 
  dplyr::group_by(conflict_id,dyad_id) %>% 
  dplyr::mutate_at(c("ext_sup","ext_sup_gov","ext_sup_rebel"), 
                   list(lag1 = ~dplyr::lag(.,n=1L))) %>% 
  ungroup()

summary(paperdata)

## DV: 
# ext_count
#ext_x_count
#ext_p
# ext_y
#ext_y_coun


dv <- "ext_sup"
## controls
ctrs <- c("polity2_lag1","gdp_log_lag1", "military_log_lag1", "population_log_lag1",
          "bd_best_lag1", "deaths_civilians_lag1", "factor(gwno_loc)", "factor(year)")
ivs_1 <- c("troops_backward_lag1")
ivs_2 <- c("police_backward_lag1")
ivs_3 <- c("observers_backward_lag1")
ivs_4 <- c("totals_backward_lag1")


f_1 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_1, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_2 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_2, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_3 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_3, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_4 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_4, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
ivs <- c("troops_backward_lag1","police_backward_lag1","observers_backward_lag1","totals_backward_lag1")


## modeling
model_ext_sup <- list()
for (i in 1:length(ivs)){
  model_ext_sup[[i]] <- glm(eval(parse(text = paste0('f_', i))), 
                            data = paperdata, family = binomial(link = "logit"))
}







p1 <- marginaleffect_logit(ModelResults = model_ext_sup, n.sim = 1000, 
                           varname = c("troops_backward_lag1","police_backward_lag1",
                                       "observers_backward_lag1","totals_backward_lag1"),
                           data = paperdata,
                           val1 = 0, val2 =1, clusterid = "gwno_loc") 

p1 + scale_x_continuous(breaks = 1:4,
                        labels = parse(text =c( 
                          "维和部队",
                          "维和警察",
                          "维和观察员",
                          "所有类型")))+ 
  labs(title="因变量：是否存在外部支持",
       tag = "图2a ")+
  theme(plot.tag.position = "bottom")
# ggsave("Results/figures/fd_ext_sup.pdf", width = 6, height = 4)
ggsave("Replication/figures/figure2-a.bmp", width = 6, height = 4)
ggsave("Replication/figures/figure2-a.pdf", width = 6, height = 4)



### 支持政府或反政府


dv <- "ext_sup_gov"
## controls
ctrs <- c(#"ext_sup_gov_lag1",
  "polity2_lag1","gdp_log_lag1", "military_log_lag1", "population_log_lag1",
  "bd_best_lag1", "deaths_civilians_lag1", "factor(gwno_loc)", "factor(year)")
ivs_1 <- c("troops_backward_lag1")
ivs_2 <- c("police_backward_lag1")
ivs_3 <- c("observers_backward_lag1")
ivs_4 <- c("totals_backward_lag1")


f_1 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_1, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_2 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_2, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_3 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_3, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_4 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_4, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
ivs <- c("troops_backward_lag1","police_backward_lag1","observers_backward_lag1","totals_backward_lag1")


## modeling
model_ext_supgov <- list()
for (i in 1:length(ivs)){
  model_ext_supgov[[i]] <- glm(eval(parse(text = paste0('f_', i))), 
                               data = paperdata, family = binomial(link = "logit"))
}


dv <- "ext_sup_rebel"
## controls
ctrs <- c(#"ext_sup_rebel_lag1",
  "polity2_lag1","gdp_log_lag1", "military_log_lag1", "population_log_lag1",
  "bd_best_lag1", "deaths_civilians_lag1", "factor(gwno_loc)", "factor(year)")
ivs_1 <- c("troops_backward_lag1")
ivs_2 <- c("police_backward_lag1")
ivs_3 <- c("observers_backward_lag1")
ivs_4 <- c("totals_backward_lag1")


f_1 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_1, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_2 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_2, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_3 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_3, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_4 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_4, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
ivs <- c("troops_backward_lag1","police_backward_lag1","observers_backward_lag1","totals_backward_lag1")


## modeling
model_ext_suprebel <- list()
for (i in 1:length(ivs)){
  model_ext_suprebel[[i]] <- glm(eval(parse(text = paste0('f_', i))), 
                                 data = paperdata, family = binomial(link = "logit"))
}


# p1 <- marginaleffect_logit(ModelResults = model_ext_suprebel, n.sim = 1000, 
#                            varname = c("troops_backward_lag1","police_backward_lag1",
#                                        "observers_backward_lag1","totals_backward_lag1"),
#                            data = paperdata,
#                            val1 = 0, val2 =1, clusterid = "gwno_loc") 
# 
# p1 + scale_x_continuous(breaks = 1:4,
#                         labels = parse(text =c( 
#                           "维和部队",
#                           "维和警察",
#                           "维和观察员",
#                           "所有类型")))+ ggtitle("是否存在外部支持政府")
# 
# 


p1_1 <- marginaleffect_logit2(ModelResults =  model_ext_supgov , n.sim = 1000, 
                              varname = c("troops_backward_lag1","police_backward_lag1",
                                          "observers_backward_lag1","totals_backward_lag1"),
                              data = paperdata,
                              dv = "ext_sup_gov",
                              val1 = 0, val2 =1, clusterid = "gwno_loc") 


p1_2 <- marginaleffect_logit2(ModelResults =  model_ext_suprebel, n.sim = 1000, 
                              varname = c("troops_backward_lag1","police_backward_lag1",
                                          "observers_backward_lag1","totals_backward_lag1"),
                              data = paperdata,
                              dv = "ext_sup_rebel",
                              val1 = 0, val2 =1, clusterid = "gwno_loc") 


df_ext_sup <- bind_rows(p1_1, p1_2)



p = ggplot(df_ext_sup,aes(color = dv)) + 
  geom_pointrange(aes(x = id, y =median, ymin =low,
                      ymax =high, shape = dv),
                  fatten = 7,
                  lwd = 1.5, position = position_dodge(width = 1/2),
                  fill="white")+
  geom_hline(yintercept = 0, lty = 2) +
  geom_linerange(aes(x = id, ymin = low,
                     ymax = high),
                 lwd = 1.5, position = position_dodge(width = 1/2)) +

  geom_text(aes(x = id +0.2 ,
                y = median,
                label = format(round(median, 2))),position = position_dodge(width = 1/2)) +
  theme_bw() + xlab("") + ylab("预测概率变化（边际效应）") + 
  theme(legend.position="bottom",
        legend.title=element_blank(),
        legend.text = element_text(colour="black", size=16),
        axis.text= element_text(size=16),
        text = element_text(size=16, family = 'KaiTi',face = "bold"),
        plot.title = element_text(hjust = .5, size = 16))+
    labs(title="因变量：是否存在外部支持",
      tag = "图2b ")+
  theme(plot.tag.position = "bottom")+
  scale_colour_manual(values = c("black", "black"), labels = c("支持政府","支持反叛组织"))+
  scale_shape_manual(values = c(15,23), labels = c("支持政府","支持反叛组织"))+
  guides(color = guide_legend(nrow = 1))+
  scale_x_continuous(breaks = 1:4,
                     labels = parse(text =c( 
                       "维和部队",
                       "维和警察",
                       "维和观察员",
                       "所有类型")))
#ggsave("Results/figures/fd_ext_sup_type.pdf", plot = p, width = 6, height = 4)
ggsave("Replication/figures/figure2-b.bmp", width = 6, height = 4)
ggsave("Replication/figures/figure2-b.pdf", width = 6, height = 4)




save(model_ext_sup,model_ext_supgov,model_ext_suprebel, file = "Results/Figure1nodels.RData")

modSD = cluster_se(ModelResults = c(model_ext_sup,model_ext_supgov,model_ext_suprebel), data = paperdata, 
                   clusterid = "gwno_loc") 

modPvalue = cluster_pvalue(ModelResults = c(model_ext_sup,model_ext_supgov,model_ext_suprebel), data = paperdata, 
                           clusterid = "gwno_loc") 

screenreg(c(model_ext_sup,model_ext_supgov,model_ext_suprebel ),override.se =  modSD,override.pvalues = modPvalue )
## table-1
htmlreg(c(model_ext_sup,model_ext_supgov,model_ext_suprebel ), file = "Replication/tables/ext_sup.doc",
        stars = c(0.01, 0.05, 0.1),
        override.se =  modSD,
        override.pvalues = modPvalue,
        caption = "Logit回归结果：维和行动对外部支持的影响",
        caption.above = TRUE,
        label = "tab:tab5",
        scalebox = 0.7,
        custom.header = list("外部支持" = 1:4,
                             "支持政府" = 5:8,
                             "支持反叛组织" = 9:12),
        # custom.model.names = c("M1:troop", "M2:police", "M3:observers", "M4:any",
        #                        "M5:troop", "M6:police", "M7:observers", "M8:any",
        #                        "M9:troop", "M10:police", "M11:observers", "M12:any"),
        custom.coef.map = list("troops_backward_lag1" = "维和部队",
                               "police_backward_lag1" = "维和警察",
                               "observers_backward_lag1" ="维和观察员",
                               "totals_backward_lag1" = "所有类型",
                               "polity2_lag1" = "政体类型",
                               "gdp_log_lag1" = "国民生产总值（对数）",
                               "military_log_lag1" = "军队人数（对数）",
                               "population_log_lag1" = "人口总数(对数)",
                               "bd_best_lag1"= "战斗相关死亡人数",
                               "deaths_civilians_lag1" = "平民死亡人数"),
        custom.gof.rows = list("国家固定效应" = c("是", "是", "是", "是",
                                            "是", "是", "是", "是",
                                            "是", "是", "是", "是"),
                               "年份固定效应" = c("是", "是", "是", "是",
                                            "是", "是", "是", "是",
                                            "是", "是", "是", "是")),
        digits = 3)


#############################################################